﻿/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Jurgen Van Gael
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */

using System;

namespace dnAnalytics.Statistics.Distributions
{
    /// <summary>
    /// The multinomial distribution is a distribution over arrays of integers. The distribution is parameterized by
    /// a vector of probabilities. In dnAnalytics, the multinomial distribution is parameterized by an array of
    /// nonnegative numbers. This array does not need to be normalized (as this is often impossible using floating
    /// point arithmetic).
    /// </summary>
    /// <remarks>The distribution will use the <see cref="System.Random"/> by default. 
    /// Users can set the random number generator by using the <see cref="RandomNumberGenerator"/> property.
    /// 
    /// The statistics classes will check all the incoming parameters whether they are in the allowed
    /// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters
    /// to false, all parameter checks can be turned off.</remarks>
    public class Multinomial
    {
        // Stores the normalized multinomial probabilities.
        private readonly double[] mP;

        /// <summary>
        /// Constructs a multinomial distribution.
        /// </summary>
        /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized 
        /// as this is often impossible using floating point arithmetic.</param>
        /// <exception cref="ArgumentException">If any of the probabilities are negative or do not sum to one.</exception>
        public Multinomial(double[] p)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(p);
            }

            mP = p;
            RandomNumberGenerator = new System.Random();
        }

        /// <summary>
        /// Generate a multinomial distribution from histogram <paramref name="h"/>. The distribution will
        /// not be automatically updated when the histogram changes.
        /// </summary>
        public Multinomial(Histogram h)
        {
            // The probability distribution vector.
            mP = new double[h.BinCount];

            // Fill in the distribution vector.
            for (int i = 0; i < h.BinCount; i++)
            {
                mP[i] = h[i];
            }

            RandomNumberGenerator = new System.Random();
        }

        /// <summary>
        /// Check the parameters of the distribution.
        /// </summary>
        /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized 
        /// as this is often impossible using floating point arithmetic.</param>
        /// <exception cref="ArgumentException">If any of the probabilities are negative or do not sum to one.</exception>
        public static void CheckParameters(double[] p)
        {
            for (int i = 0; i < p.Length; i++)
            {
                if (p[i] < 0.0)
                {
                    throw new ArgumentException("The entries of p are probabilities and hence cannot be 0.");
                }
            }
        }

        /// <summary>
        /// Returns the parameters of the multinomial distribution: this can be either an array of ratios or
        /// the actual normalized probabilities.
        /// </summary>
        /// <param name="normalized">Whether to normalize the parameters to (within floating point constraints) 1 or not.</param>
        /// <returns></returns>
        public double[] P(bool normalized)
        {
            if (!normalized)
            {
                return mP;
            }

            double sum = 0.0;
            for (int i = 0; i < mP.Length; i++)
            {
                sum += mP[i];
            }
            double[] p = (double[]) mP.Clone();
            for (int i = 0; i < mP.Length; i++)
            {
                p[i] /= sum;
            }
            return p;
        }

        /// <summary>
        /// Samples one multinomial distributed random variable; sometimes known as the Multinoulli distribution.
        /// </summary>
        public int Sample()
        {
            return Sample(RandomNumberGenerator, mP);
        }

        /// <summary>
        /// Samples an array of Bernoulli distributed random variables.
        /// </summary>
        /// <param name="n">The number of variables needed.</param>
        public int[] Sample(int n)
        {
            return Sample(RandomNumberGenerator, n, mP);
        }

        /// <summary>
        /// Samples one multinomial distributed random variable; sometimes known as the Multinoulli distribution.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized 
        /// as this is often impossible using floating point arithmetic.</param>
        public static int Sample(System.Random rnd, double[] p)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(p);
            }

            // Compute the cumulative density of p.
            double[] cp = (double[]) p.Clone();
            for (int i = 1; i < p.Length; i++)
            {
                cp[i] += cp[i - 1];
            }

            double u = rnd.NextDouble()*cp[cp.Length - 1];
            int idx = 0;
            while (u > cp[idx])
            {
                idx++;
            }
            return idx;
        }

        /// <summary>
        /// Samples an array of Bernoulli distributed random variables.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="n">The number of variables needed.</param>
        /// <param name="p">An array of nonnegative ratios: this array does not need to be normalized 
        /// as this is often impossible using floating point arithmetic.</param>
        public static int[] Sample(System.Random rnd, int n, double[] p)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(p);
            }

            // Compute the cumulative density of p.
            double[] cp = (double[]) p.Clone();
            for (int i = 1; i < p.Length; i++)
            {
                cp[i] += cp[i - 1];
            }

            int[] arr = new int[n];
            for (int i = 0; i < n; i++)
            {
                double u = rnd.NextDouble()*cp[cp.Length - 1];
                int idx = 0;
                while (u > cp[idx])
                {
                    idx++;
                }
                arr[i] = idx;
            }

            return arr;
        }

        /// <summary>
        /// Gets or sets the random number generator.
        /// </summary>
        /// <value>The random number generator used to generate a random sample.</value>
        public System.Random RandomNumberGenerator { get; set; }
    }
}